Citation: Elio J et al. Application of airborne radiometric surveys for large-scale geogenic radon potential classification. Journal of the European Radon Association 2020, 1: 4358, http://dx.doi.org/10.35815/radon.v1.4358.

Corresponding author: Javier Elío; Centre for the Environment, Trinity College, Dublin 2, Ireland. ORCID: 0000-0003-0624-2345; - .

Abstract:

Background: Indoor radon represents an important health issue to the general population. Therefore, accurate radon risk maps help public authorities to prioritise areas where mitigation actions should be implemented. As the main source of indoor radon is the soil where the building is constructed, maps derived from geogenic factors ([e.g. geogenic radon potential [GRP]) are viewed as valuable tools for radon mapping.

Objectives: A novel indirect method for estimating the GRP at national/regional level is presented and evaluated in this article.

Design: We calculate the radon risk solely based on the radon concentration in the soil and on the subsoil permeability. The soil gas radon concentration was estimated using airborne gamma-ray spectrometry (i.e. equivalent uranium [eU]), assuming a secular equilibrium between eU and radium (226Ra). The subsoil permeability was estimated based on groundwater subsoil permeability and superficial geology (i.e. quaternary geology) by assigning a permeability category to each soil type (i.e. low, moderate or high). Soil gas predictions were compared with in situ radon measurements for representative areas, and the resulting GRP map was validated with independent indoor radon data.

Results: There was good agreement between soil gas radon predictions and in situ measurements, and the resultant GRP map identifies potential radon risk areas. Our model shows that the probability of having an indoor radon concentration higher than the Irish reference level (200 Bq m-3) increases from c. 6% (5.2% – 7.1%) for an area classified as low risk, to c. 9.7% (9.1% – 10.5%) for moderate-low risk areas, c. 14% (13.4% – 15.3%) for moderate-high risk areas and c. 26% (24.5% – 28.6%) for high-risk areas.

Conclusions: The method proposed here is a potential alternative approach for radon mapping when airborne radiometric data (i.e. eU) are available.

Data

Ireland (RoI and NI)

Read grid cells of 1km x 1 km (same grid cells where Ireland has population data just in case we would need to link the RP radon map with population density).

  Grids1km <- read_sf("Rdata/Grids1km/IR_NI_Grids_1km.shp")
    
  Ireland <- Grids1km %>%
    group_by(Id) %>%
    summarise(N = n()) %>%
    ungroup()

Airborne (eU)

We have aggregated the initial data (with a resolution of 50 x 50 m; i.e. Tellus radiometric merged grd, tiff and gxf (981.9 MB) - Tellus) in grid cells of 1x1 km (after running a separate script, i.e. Average_eU_grid_cells_1km2). Here we read directly the aggregated map (eU_Grids1km), and we use the initial map (eU_50m) only for comparison (Figure 3).

  eU<-"Rdata/Radiometrics/RAD_MERGED_GRIDS/GXF/Tellus_RAD_eU_MERGE_2019.gxf"
  eU_50m  <- read_stars(eU) %>% st_set_crs(st_crs(Grids1km)) 
     
  eU_Grids1km <- read_sf("Rresults/eU_Grids1km.shp")

  eU_50m_tbl <- eU_50m %>% 
  #as.tbl_cube() %>%
  as.data.frame() %>%
  drop_na()  
  
  cbind(tibble(Resolution = c("50 m","1 km")),
        rbind(summary(eU_50m_tbl$Tellus_RAD_eU_MERGE_2019.gxf),
              summary(eU_Grids1km$eU_AM)
        )
  ) %>%
    mutate_if(is.numeric, ~round(., 2)) %>%
    kbl() %>%
    kable_styling()
Resolution Min. 1st Qu. Median Mean 3rd Qu. Max.
50 m 0 0.33 0.72 0.71 0.99 10.0
1 km 0 0.39 0.73 0.73 0.97 6.6

Radium predictions (Ra-226)

  • Radium (Ra-226) activity in Bq/kg
  • Equilibrium with eU (1 ppm eU = 12.35 Bq/kg)
  eU_Grids1km <- eU_Grids1km %>%
    mutate(AM_Ra226 = 12.35*eU_AM,
           SD_Ra226 = 12.35*eU_SD) 
  ggplot() +
    geom_sf(data = Ireland) +
    geom_sf(data = eU_Grids1km, aes(fill = AM_Ra226), colour = NA) +
    scale_fill_gradient2(
    name = expression(""^226*Ra * " [ " * "Bq "* kg^-1 * " ]"),
    midpoint = 30, 
    low   = "#0000AA",
    mid   = "#AAFFAA",
    high  = "#AA0000") +
    ggtitle("Airborne radiometrics (Grids 1km x 1km)") 

All-Ireland Quaternary map (scale 1:500.000)

  QG <- read_sf("Rdata/GSI_Quaternary_500k/QG_500k.shp") %>% 
    mutate(Sediment_5 = as_factor(Sediment_5)) %>%
    group_by(Sediment_5) %>%
    summarise(N = n(),
              area_sq_km = sum(area_sq_km),
              label = first(label)) %>%
    ungroup()

Aquifer Type

AQ <- read_sf("Rdata/Groundwater/IRL_AQUIFER_BEDROCK_ITM.shp") %>%
    rmapshaper::ms_simplify(keep = 0.01, keep_shapes = TRUE) %>%
    st_make_valid()
  
  Karst <- AQ %>% 
    group_by(AQUIFERCAT) %>%
    summarize(AQUIFER_DE = first(AQUIFER_DE)) %>%
    ungroup() %>%
    filter(AQUIFERCAT != "Lake" & AQUIFERCAT != "Unclas") %>%
    mutate(Karst = case_when(
      AQUIFERCAT == "Lk"     ~ "TRUE",
      AQUIFERCAT == "LkNI"   ~ "TRUE",
      AQUIFERCAT == "Ll"     ~ "FALSE",                         
      AQUIFERCAT == "LlNI"   ~ "FALSE",                 
      AQUIFERCAT == "Lm"     ~ "FALSE",                
      AQUIFERCAT == "LmNI"   ~ "FALSE",              
      AQUIFERCAT == "NI_Un*" ~ "FALSE",                  
      AQUIFERCAT == "PINIl"  ~ "FALSE",       
      AQUIFERCAT == "Pl"     ~ "FALSE",      
      AQUIFERCAT == "PlNI"   ~ "FALSE",     
      AQUIFERCAT == "Pu"     ~ "FALSE",    
      AQUIFERCAT == "PuNI"   ~ "FALSE",   
      AQUIFERCAT == "Rf"     ~ "FALSE",  
      AQUIFERCAT == "Rf/Rk"  ~ "TRUE" , 
      AQUIFERCAT == "RfNI"   ~ "FALSE",
      AQUIFERCAT == "Rk"     ~ "TRUE",
      AQUIFERCAT == "Rkc"    ~ "TRUE",
      AQUIFERCAT == "RkcNI"  ~ "TRUE",
      AQUIFERCAT == "Rkd"    ~ "TRUE",
      AQUIFERCAT == "RkNI"   ~ "TRUE",
      TRUE ~ as.character(AQUIFERCAT))
    ) %>%
    group_by(Karst) %>%
    summarise(N = n()) %>%
    ungroup()

Subsoil permeability

Load SP_GW_naQG.shp file; output of the R-script: SP_Groundwater_recharge.R.

  SP <- read_sf("Rresults/Soil_Permeability/SP_GW_naQG.shp") %>%
    rename(Perm = PERM) %>%
    mutate(Perm = factor(Perm,
                         levels = c("H", "M", "L"),
                         ordered = TRUE)) 

Rasterize: convert Subsoil Permeability shapefile (SP) to a raster file with a resolution of 1km x 1 km).

  1. Make grids equal to the cell grids of eU data (1km x 1km)
  grd_res_x <- 1000 # m
  grd_res_y <- 1000 # m
  grd <- st_as_stars(st_bbox(eU_Grids1km),
     nx = (st_bbox(eU_Grids1km)[3] - st_bbox(eU_Grids1km)[1]) / grd_res_x,
     ny = (st_bbox(eU_Grids1km)[4] - st_bbox(eU_Grids1km)[2]) / grd_res_y,
     xlim = c(st_bbox(eU_Grids1km)[1], st_bbox(eU_Grids1km)[3]),
     ylim = c(st_bbox(eU_Grids1km)[2], st_bbox(eU_Grids1km)[4]),
     values = NA_real_)
  1. Rasterize
SP_Raster <- stars::st_rasterize(SP["Perm"], grd) %>%
  rename(Perm_class = Perm) %>% 
  mutate(Perm_class = factor(Perm_class, 
                             levels = c("H", "M", "L"),
                             ordered = TRUE)) %>%
  st_crop(eU_Grids1km)

Results

Soil-gas radon predictions based on airborne radiometrics (eU)

  1. Soil Parameters:
  • e ~ N(0.29, 0.03); Yu et al 1993: Range approx. from 0.13 to 0.37
  • d ~ N(1.35, 0.06); Yu et al 1993: Range dry density 1.1 - 1.6; Creamer et al 2016 (Irish soils) bulk density mainly 1.1-1.3 and >1.4 g m-3
  • n ~ N(0.30, 0.07); Yu et al 1993: Effective Porosity range approx. from 0.02 to 0.53
  • Water influence (SF) ~ uniform distribution on the interval from min (0.4) to max (0.6)
  set.seed(1) # Make the analysis reproducible
  
  e  <- rnorm(1000, mean = 0.29, sd = 0.03)
  d  <- rnorm(1000, mean = 1.35, sd = 0.06)
  n  <- rnorm(1000, mean = 0.30, sd = 0.07)
  SF <- runif(1000, min = 0.4, max = 0.6)
  id <- "g90164" 
  Ra <- rnorm(1000,
              mean = eU_Grids1km[eU_Grids1km$SP_ID == id, ]$AM_Ra226,
              sd   = eU_Grids1km[eU_Grids1km$SP_ID == id, ]$SD_Ra226
              )
  
  par(mfrow = c(2,3))
  hist(e, col = 2, breaks = 15, main = "Emanation coefficient")
  hist(d, col = 2, breaks = 15, main = "Soil density",
       xlab = expression(d * " [" *g *" "* cm^-3 * "]"))
  hist(n, col = 2, breaks = 15, main = "Effective porosity")
  hist(SF, col = 2, breaks = 15, main = "Water saturation")
  hist(Ra, col = 2, breaks = 15, main = "Ra-226")

  1. Simulations

Radon concentration in soil-gas:

\[C_{Rn} [kBq m^{-3}] = \frac{C_{Ra} [Bq kg^{-1}] \times e \times \rho [g m^{-3}]}{n} \times \frac{1}{1 - S_F + S_FK_{W/Air}}\]

  set.seed(1) # Make the analysis reproducible

  rep <- 1000
  SG_MCS <- matrix(NA, nrow = dim(eU_Grids1km)[1], ncol = rep+1)
  for (i in 1:rep) {
    Ra <- rnorm(dim(eU_Grids1km)[1],
                eU_Grids1km$AM_Ra226,
                eU_Grids1km$SD_Ra226)
    e  <- rnorm(dim(eU_Grids1km)[1], 0.29, 0.03)
    d  <- rnorm(dim(eU_Grids1km)[1], 1.35, 0.06)
    n  <- rnorm(dim(eU_Grids1km)[1], 0.30, 0.07)
    SF <- runif(dim(eU_Grids1km)[1], 0.40, 0.60)
    SG_MCS[ ,i+1] <- (Ra*e*d/n)*(1/(1-SF+(0.25*SF)))
  }
  SG_MCS <- as.data.frame(SG_MCS)
  names(SG_MCS) <- c("SP_ID", c(1:rep))
  SG_MCS$SP_ID <- eU_Grids1km$SP_ID
   
  SG_MCS <- as_tibble(SG_MCS)
  SG_MCS
  1. Percentiles:
  SG_MCS_T <- t(SG_MCS[,2:rep])
  
  SG_MCS_P95 <- matrix(NA, nrow = dim(eU_Grids1km)[1], ncol = 8)  
  for (i in 1:dim(eU_Grids1km)[1]) {
    SG_MCS_P95[i,2] <- quantile(SG_MCS_T[,i], 0.025)
    SG_MCS_P95[i,3] <- quantile(SG_MCS_T[,i], 0.05)
    SG_MCS_P95[i,4] <- quantile(SG_MCS_T[,i], 0.25)
    SG_MCS_P95[i,5] <- quantile(SG_MCS_T[,i], 0.50)
    SG_MCS_P95[i,6] <- quantile(SG_MCS_T[,i], 0.75)
    SG_MCS_P95[i,7] <- quantile(SG_MCS_T[,i], 0.95)
    SG_MCS_P95[i,8] <- quantile(SG_MCS_T[,i], 0.975)
  }
  
  SG_MCS_P95 <- as.data.frame(SG_MCS_P95)
  names(SG_MCS_P95) <- c("SP_ID","Rn_P025","Rn_P05","Rn_P25",
                         "Rn_P50","Rn_P75","Rn_P95","Rn_P975")
  SG_MCS_P95$SP_ID <- SG_MCS$SP_ID
  SG_MCS_P95 <- as_tibble(SG_MCS_P95)
  
  eU_Grids1km <- left_join(eU_Grids1km, SG_MCS_P95, by = "SP_ID")
  eU_Grids1km
  1. Soil gas radon classification based on eU (Rn predictions - percentile 75%):
  • Extremely High: >= 100
  • Very High: 70 - 100
  • High: 50 - 70
  • Moderate: 30 - 50
  • Low: 10 - 30
  • Very Low: < 10
  eU_Grids1km <- eU_Grids1km %>%
    mutate(Rn_class = case_when(
      Rn_P75 >= 100              ~ "EH",
      Rn_P75 >= 70 & Rn_P75 < 100 ~ "VH",
      Rn_P75 >= 50 & Rn_P75 <  70 ~  "H",
      Rn_P75 >= 30 & Rn_P75 <  50 ~  "M",                 
      Rn_P75 >= 10 & Rn_P75 <  30 ~  "L",                
      Rn_P75 <  10                ~ "VL",
      TRUE ~ as.character(Rn_P75)
    )
    ) %>%
    mutate(Rn_class = factor(Rn_class, 
                             levels = c("EH", "VH", "H", "M", "L", "VL"),
                             ordered = TRUE
    )
    )

Radon Potential map

  1. Merge Rn_class and SP in one dataset (by grids of 1km x 1km):
  • Create a dataframe with Rn values (Rn_class)
  • Transform SP_Raster to points (Perm value at the centroid)
  • Add points to the Grids1km (st_join - inner join)
  1. Assign values of Rn and Perm
  2. Estimate RP

\[ RP = \frac{C_Rn}{(-log_{10}(k) - 10)}\]

  1. RP classification (RP_class)
  • High: RP ≥ 30
  • Moderate-High: 22.5 <= RP < 30
  • Moderate-Low: 10 <= RP < 22.5
  • Low: RP < 10
  RP_Grids <- eU_Grids1km %>%
    select(Rn_class)
  
  SP_Points <- st_as_sf(SP_Raster, as_points = TRUE, merge = FALSE) %>%
    drop_na() %>%
    select(Perm_class) %>%
    rename(Perm = Perm_class)
  
  RP_Grids <- st_join(RP_Grids, SP_Points,
                      join = st_intersects,
                      left = FALSE
  )

  RP_Grids  <-  RP_Grids %>%
    mutate(Rn = case_when(
      Rn_class == "EH" ~ 110,
      Rn_class == "VH" ~  85,
      Rn_class == "H"  ~  65,
      Rn_class == "M"  ~  50,
      Rn_class == "L"  ~  25,
      Rn_class == "VL" ~   5
    )) %>%
    mutate(logk = case_when(
      Perm == "H" ~ 11,
      Perm == "M" ~ 12,
      Perm == "L" ~ 13
    )) %>%
    mutate(RP = (Rn/(logk - 10))) %>%
    mutate(RP_class = case_when(
      RP >= 30                 ~  "H",
      RP >= 22.5 &  RP <  30   ~ "MH",
      RP >= 10   &  RP <  22.5 ~ "ML",                 
      RP <  10                 ~  "L"
      )) %>%
    mutate(RP_class = factor(RP_class, 
                             levels = c("H", "MH", "ML", "L"),
                             ordered = TRUE)
    )

    
  write_sf(RP_Grids, "Rresults/RP_Grids1km.shp")

In-situ soil-gas measurements

  Soil_Rn <- readxl::read_xlsx("Rdata/data.xlsx",
                               sheet = "DATA (MeanReplica_NoNA)") %>%
    rename(SoilRn = "Rn-222") %>%
    select(Grid, Point, Day, X_ITM, Y_ITM, SoilRn) %>%
    st_as_sf(coords = c("X_ITM", "Y_ITM")) %>%
    st_set_crs(st_crs(eU_Grids1km)) %>%
    st_intersection(eU_Grids1km)

Summary

  SG_Sum <- Soil_Rn %>% 
    group_by(Grid) %>%
    summarize(N = n(),
              Rn_Min = min(SoilRn),
              Rn_GM = exp(mean(log(SoilRn))),
              Rn_AM = mean(SoilRn),
              Rn_Max = max(SoilRn),
              Rn_RMD = 100 * mean(abs(SoilRn - mean(SoilRn)))/mean(SoilRn),
              Ra_AM = mean(AM_Ra226),
              Ra_SD = mean(SD_Ra226),
              Rn_P025 = mean(Rn_P025),
              Rn_P50 = mean(Rn_P50),
              Rn_P75 = mean(Rn_P75),
              Rn_P975 = mean(Rn_P975)
              ) %>%
    ungroup()

  SG_Sum %>%
    mutate_if(is.numeric, ~round(., 0)) %>%
    as_tibble() %>%
    select(-geometry) %>% 
    kbl() %>%
    kable_styling()
Grid N Rn_Min Rn_GM Rn_AM Rn_Max Rn_RMD Ra_AM Ra_SD Rn_P025 Rn_P50 Rn_P75 Rn_P975
G56114 9 16 67 80 125 43 6 3 0 12 17 30
G57626 10 49 108 113 170 23 14 1 17 29 35 58
G64863 11 85 145 155 310 30 19 2 23 39 48 77
G69054 9 20 37 40 69 33 14 1 18 29 35 52
G70565 9 13 33 37 65 38 9 1 12 19 23 38
G70947 9 13 36 39 59 32 13 2 16 27 33 54
G72081 9 17 51 60 114 46 12 1 15 24 30 50
G73606 8 76 179 196 264 33 27 4 30 55 69 111
G78164 10 79 122 128 191 29 19 2 23 39 48 76
G78544 10 42 102 120 256 50 23 3 29 47 58 89
G90163 9 31 111 138 335 51 40 4 50 84 103 164
G98159 10 16 94 121 260 55 30 4 37 63 79 129
G98929 10 60 142 163 280 43 36 4 45 77 93 153

Indoor radon concentration

  InRn <- readxl::read_xls("Rdata/EPA/XY_InRn_EPA_Old_New.xls") %>%
    st_as_sf(coords = c("X_ITM", "Y_ITM")) %>%
    st_set_crs(st_crs(RP_Grids)) %>% #IRENET95/Irish Transverse Mercator
    mutate(Case = if_else(Rn > 200, TRUE, FALSE)) %>%
    rename(InRn   = Rn)
  Sum_New <- tibble(N = length(InRn[InRn$Survey=="New",]$InRn),
                  Prop = 100*mean(InRn[InRn$Survey=="New",]$Case),
                  Min = min(InRn[InRn$Survey=="New",]$InRn),
                  Q1 = quantile(InRn[InRn$Survey=="New",]$InRn, prob = 0.25),
                  Median = median(InRn[InRn$Survey=="New",]$InRn),
                  Mean = mean(InRn[InRn$Survey=="New",]$InRn),
                  Q3 = quantile(InRn[InRn$Survey=="New",]$InRn, prob = 0.75),
                  Max = max(InRn[InRn$Survey=="New",]$InRn),
                  SD = sd(InRn[InRn$Survey=="New",]$InRn),
                  GM = exp(mean(log(InRn[InRn$Survey=="New",]$InRn))),
                  GSD = exp(sd(log(InRn[InRn$Survey=="New",]$InRn)))
  )
  
  Sum_Old <- tibble(N = length(InRn[InRn$Survey=="Old",]$InRn),
                  Prop = 100*mean(InRn[InRn$Survey=="Old",]$Case),
                  Min = min(InRn[InRn$Survey=="Old",]$InRn),
                  Q1 = quantile(InRn[InRn$Survey=="Old",]$InRn, prob = 0.25),
                  Median = median(InRn[InRn$Survey=="Old",]$InRn),
                  Mean = mean(InRn[InRn$Survey=="Old",]$InRn),
                  Q3 = quantile(InRn[InRn$Survey=="Old",]$InRn, prob = 0.75),
                  Max = max(InRn[InRn$Survey=="Old",]$InRn),
                  SD = sd(InRn[InRn$Survey=="Old",]$InRn),
                  GM = exp(mean(log(InRn[InRn$Survey=="Old",]$InRn))),
                  GSD = exp(sd(log(InRn[InRn$Survey=="Old",]$InRn)))
  )
  
  Sum <- rbind(Sum_New, Sum_Old)
  
  Sum %>%
    mutate_if(is.numeric, ~round(., 2)) %>%
    kbl() %>%
    kable_styling()
N Prop Min Q1 Median Mean Q3 Max SD GM GSD
6859 11.90 5.0 31.00 55 109.71 110.00 6240 204.47 62.68 2.58
31910 12.46 5.7 33.19 59 115.32 115.68 9714 219.72 65.76 2.59

t-test

  t.test(InRn[InRn$Survey == "New",]$InRn,
         InRn[InRn$Survey == "Old",]$InRn)

  t.test(log10(InRn[InRn$Survey == "New",]$InRn),
         log10(InRn[InRn$Survey == "Old",]$InRn))

Binomial

  xtabs(~ Case , data = InRn[InRn$Survey == "New", ])
    binom.test(c(816,6043))           # Prob 11.89% (11.14% - 12.69%)
    
  xtabs(~ Case , data = InRn[InRn$Survey=="Old",])
    binom.test(c(3975,27935))           # Prob 12.45% (12.10% - 12.82%)
    
  prop.test(c(816,3975),c(816+6043,3975+27935))
  fisher.test(rbind(c(816,6043), c(3975,27935)), alternative="less")

Indoor radon vs. RP map

  InRn_RP <- st_join(InRn, RP_Grids, join = st_intersects) %>%
    drop_na() %>%
    mutate(RP_class = factor(RP_class, ordered = TRUE, 
                             levels = c("L", "ML", "MH", "H")))

# Intersection with 30m tolerance
  # InRn_RP <- st_join(InRn, RP_Grids, join = st_nn, maxdist = 30) %>%
  #   drop_na() %>%
  #   mutate(RP_class = factor(RP_class, ordered = TRUE, 
  #                            levels = c("L", "ML", "MH", "H")))
# GM and GSD for each RP classification
  SM <- InRn_RP %>%
    group_by(RP_class) %>%
    summarize(GM  = exp(mean(log(InRn))),
              GSD = exp(sd(log(InRn))) 
              ) %>%
    ungroup() %>%
    as_tibble() %>%
    mutate(RP_class = factor(RP_class, ordered = TRUE, 
                             levels = c("L", "ML", "MH", "H"))) %>%
    arrange(RP_class) %>%
    select(RP_class, GM, GSD)
  
# Cases (InRn > 200 Bq/m3)
  CM <- xtabs(~ Case + RP_class , data = InRn_RP) %>%
    as_tibble() %>%
    spread(key = Case, value = n) %>%
    mutate(RP_class = factor(RP_class, ordered = TRUE, 
                             levels = c("L", "ML", "MH", "H"))) %>%
    arrange(RP_class) %>%
    rename(No_Case = "FALSE", 
           Case =  "TRUE") %>%
    mutate(Total = No_Case + Case)

# Binomial distribution
  BTM <- matrix(NA, 4, 3)
    for (i in 1:4) {
      bt <- binom.test(CM[i,]$Case, CM[i,]$Total)
      BTM[i, 1] <- 100 * bt$estimate
      BTM[i, 2] <- 100 * bt$conf.int[1]
      BTM[i, 3] <- 100 * bt$conf.int[2]
    }
  BTM <- BTM %>%
    as_data_frame() %>%
    rename(Prob = V1,
           LCI  = V2,
           UCI  = V3)
   
# Total
  TM <- tibble(RP_class = "Total",
               GM       = exp(mean(log(InRn_RP$InRn))),
               GSD      = exp(sd(log(InRn_RP$InRn))),
               No_Case  = length(filter(InRn_RP, Case == FALSE)$Case),
               Case     = length(filter(InRn_RP, Case == TRUE)$Case),
               Total    = length(filter(InRn_RP)$Case)
  ) %>%
    mutate(Prob = 100 * binom.test(Case, Total)$estimate,
           LCI  = 100 * binom.test(Case, Total)$conf.int[1],
           UCI  = 100 * binom.test(Case, Total)$conf.int[2])

Table_6 <- cbind(SM, select(CM, -RP_class), BTM) %>%
    rbind(TM) %>%
    mutate_if(is.numeric, ~round(., 2)) 

Table_6 %>%
    kbl() %>%
    kable_styling()
RP_class GM GSD No_Case Case Total Prob LCI UCI
L 47.50 2.36 2191 142 2333 6.09 5.15 7.13
ML 57.90 2.48 6620 716 7336 9.76 9.09 10.46
MH 70.68 2.65 4527 757 5284 14.33 13.39 15.30
H 98.65 3.01 1349 486 1835 26.49 24.48 28.57
Total 63.58 2.63 14687 2101 16788 12.51 12.02 13.02

Figures

Figure 1

  1. Example of the predicted soil-gas radon concentration in the grid G90164; b) percentiles 2.5%, 50%, 75%, and 97.5% of the simulated values for all grids covered by the Tellus gamma-ray spectrometry airborne radiometrics
Row <- as.numeric(rownames(SG_MCS_P95[SG_MCS_P95$SP_ID == id, ]))
  
cuts <- data.frame(Per = c("2.5th", "50th", "75th", "97.5th"),
                   vals = c(SG_MCS_P95[SG_MCS_P95$SP_ID == id, ]$Rn_P025,
                            SG_MCS_P95[SG_MCS_P95$SP_ID == id, ]$Rn_P50,
                            SG_MCS_P95[SG_MCS_P95$SP_ID == id, ]$Rn_P75,
                            SG_MCS_P95[SG_MCS_P95$SP_ID == id, ]$Rn_P975
                   ),
                   stringsAsFactors = FALSE)

dat_ex <- SG_MCS[SG_MCS$SP_ID == id, ] %>%
  select(-1) %>%
  t() %>%
  as_tibble()
  
p1 <- ggplot(data = dat_ex, aes(V1)) +
  geom_histogram(fill = "gray", color = "#e9ecef") +
  xlab(expression(""^222*Rn * " [ " * "Bq "* m^-3 * " ]")) +
  geom_vline(data = cuts,
             aes(xintercept = vals, 
                 colour = Per),
             linetype = "longdash",
             size = 1, 
             show.legend = T,
             key_glyph = draw_key_smooth) +
  scale_color_manual(name = "Percentiles",
                     values = c("purple",
                                "green",
                                "blue",
                                "red")) +
  ggtitle(paste0("A) ", "predicted concentrations in the grid ",
                 SG_MCS_P95[SG_MCS_P95$SP_ID == id, ]$SP_ID)) +
  theme(legend.key = element_rect(fill = NA),
        legend.background = element_rect(color = "gray",
                                         linetype = "solid"),
        legend.position = c(.95, .95),
        legend.justification = c("right", "top"),
        legend.box.just = "right",
        legend.margin = margin(6, 6, 6, 6)) 
  

breaks <- c(1, 10, 100, 1000)
minor_breaks <- rep(1:9, 21)*(10^rep(-10:10, each = 9))
p2 <- ggplot(eU_Grids1km) +
  geom_hline(yintercept = c(10, 30, 50, 70, 100),
             linetype = "dashed", color = "darkgrey", size = 0.5) +
  geom_point(aes(x = AM_Ra226, y = Rn_P025),
             color = "purple", size = 0.1) +
  geom_point(aes(x = AM_Ra226, y = Rn_P50),
             color = "green", size = 0.1) +
  geom_point(aes(x = AM_Ra226, y = Rn_P75),
             color = "blue", size = 0.1) +
  geom_point(aes(x = AM_Ra226, y = Rn_P975),
             color = "red", size = 0.1) +
  ggtitle("B) Rn-222 soil-gas predictions (Grids 1km x 1km)") + 
  xlab(expression("Mean "*""^226*Ra * " [ " * "Bq "* kg^-1 * " ]")) + 
  ylab(expression(""^222*Rn * " [ " * "kBq "* m^-3 * " ]")) +
  scale_y_log10(breaks = breaks,
                minor_breaks = minor_breaks,
                limits = c(1,1000)) 

  ggsave(file = "Rresults/Figure_1.png",
         gridExtra::arrangeGrob(p1, p2, nrow = 1),
       width = 30,
       height = 10,
       units = "cm" )

Figure 2

Equivalent uranium concentration (ppm)

  #round(BAMMtools::getJenksBreaks(eU_Grids1km$eU_AM, 10, subset = NULL), 1)

  colpal <- colorRamps::matlab.like(9)
  #colpal <- viridis::viridis(9)
  
  grDevices::png("Rresults/Figure_2a.png",
                 width  = 550,
                 height = 650,
                 unit = "px") 

  nf <- layout(matrix(1:2, ncol = 2), widths = c(0.8, lcm(2)))
 {brks <- c(0.0, 0.3, 0.5, 0.7, 0.9, 1.2, 1.6, 2.3, 3.6, 10.2)
   plot(st_geometry(Ireland),
        reset = FALSE,
        main = "a) eU [ppm] - Resolution 50m")
   plot(eU_50m,
         main = "",
         col = colpal, 
         breaks = brks,
         at = brks,
         add = TRUE,
         reset = FALSE,
         )
   .image_scale(brks,
               col = colpal, 
               key.length = lcm(10),
               key.pos = 4,
               breaks = brks,
               at = brks, 
               add.axis = F)
  axis(4,at = brks, las = 2, cex.axis = 0.75)
  }
  dev.off()
  
  grDevices::png("Rresults/Figure_2b.png",
                 width  = 550,
                 height = 650,
                 unit = "px") 
  
  nf <- layout(matrix(1:2, ncol = 2), widths = c(0.8, lcm(2)))
  {brks <- c(0.0, 0.3, 0.5, 0.7, 0.9, 1.2, 1.6, 2.3, 3.6, 6.7)
  plot(st_geometry(Ireland),
       reset = FALSE, 
       main = "b) eU [ppm] - Resolution 1km",)
  plot(eU_Grids1km["eU_AM"], border=NA,
       main = "",
       pal = colpal, 
       breaks = brks,
       at = brks,
       add = TRUE,
       reset = FALSE,
       )
  .image_scale(brks,
               col = colpal, 
               key.length = lcm(10),
               key.pos = 4,
               breaks = brks,
               at = brks, 
               add.axis = FALSE)
  axis(4,at = brks, las = 2, cex.axis = 0.75)
  }
  dev.off()

Figure 3

Subsoil permeability map of the island of Ireland (based on the All-Ireland Quaternary map)

  Karst <- Karst %>% mutate(Bedrock = if_else(Karst == TRUE,
                                              "Bedrock (Karstified)",
                                              "Bedrock (No Karstified)"))
  
  QG_Colors <- c("Bedrock (Karstified)"    = "darkgrey",
                 "Bedrock (No Karstified)" = "lightgrey",
                 "Alluvium" = "#ffff89",
                 "Glaciofluvial/glaciolacustrine sand and gravel" = "#CCCC33",
                 "Glaciomarine sediments" = "#CC0000",
                 "Lacustrine sediments" = "yellow",
                 "Marine and estuarine deposits" = "#f4efe4",
                 "Peat" = "#ac7f50",
                 "Slope deposits" = "#996699",
                 "Till derived from acidic volcanic rocks" = "#FF9999",
                 "Till derived from basic igneous rocks" = "#ffc1b7",
                 "Till derived from Devonian sandstones" = "#CD5C5C",
                 "Till derived from granites" = "#ff3317",
                 "Till derived from limestones" = "#bddbf1",
                 "Till derived from metamorphic rocks" = "#a7a7ff",
                 "Till derived from sandstones" = "#ffcb23",
                 "Till derived from sandstones and shales" = "#888888",
                 "Wind blown sand" = "#f1e5df" 
  )
                  
  p1 <- ggplot() +
    geom_sf(data = Karst, aes(fill = Bedrock), colour = NA) +
    geom_sf(data = QG, aes(fill = Sediment_5), colour = NA) +
    scale_fill_manual(name = "Deposits", values = QG_Colors) +
    ggtitle("All Ireland Quaternary Map 500k") 

  p2 <- ggplot() +
    geom_sf(data = SP, aes(fill = Perm), colour = NA) +
    scale_fill_manual(name = "Perm", 
                      labels = c("High", "Moderate", "Low"),
                      values = c("red", "yellow3", "blue")) +
    ggtitle("Subsoil permeability based on GW, QG and AT") 

  ggsave(file = "Rresults/Figure_3.png",
         ggpubr::ggarrange(p1, p2, nrow = 1, align = "v"),
         #ridExtra::arrangeGrob(p1, p2, nrow = 1),
       width = 38,
       height = 14,
       units = "cm")

Figure 4

Indoor radon measurements from old survey (Jan. 1992 – Feb. 2013; N = 31,910) and the new survey (Feb. 2013 – Jul. 2017; N = 6,859). Grey area: part covered by the Tellus airborne survey illustrated in this study. Note that no geolocated indoor radon data are available for Northern Ireland, as these are held by Public Health England and are not available for this study.

  p1 <- ggplot() +
    geom_sf(data = Ireland) +
    geom_sf(data = RP_Grids, colour = "grey") +
    geom_sf(data = filter(InRn, Survey == "Old"),
            aes(fill = Case, color = Case),
            size = 0.1,
            alpha = 0.75) +
    scale_fill_manual(name = "Bq/m3",
                      labels = c("< 200",
                                 ">= 200"),
                      values = c("black",
                                 "red")) +
    scale_color_manual(name = "Bq/m3",
                       labels = c("< 200",
                                  ">= 200"),
                       values = c("black",
                                  "red")) +
    guides(colour = guide_legend(override.aes = list(size = 2))) +
    labs(title = "Period Jan. 1992 - Feb. 2013") 

  p2 <- ggplot() +
    geom_sf(data = Ireland) +
    geom_sf(data = RP_Grids, colour = "grey") +
    geom_sf(data = filter(InRn, Survey == "New"),
            aes(fill = Case, color = Case),
            size = 0.1,
            alpha = 0.75) +
    scale_fill_manual(name = "Bq/m3",
                      labels = c("< 200",
                                 ">= 200"),
                      values = c("black",
                                 "red")) +
    scale_color_manual(name = "Bq/m3",
                       labels = c("< 200",
                                  ">= 200"),
                       values = c("black",
                                  "red")) +
    guides(colour = guide_legend(override.aes = list(size = 2))) +
    labs(title = "Period Feb. 2013 - Jul. 2017")  

  ggsave(file = "Rresults/Figure_4.png",
         gridExtra::arrangeGrob(p1, p2, nrow = 1),
         width = 28,
         height = 14,
         units = "cm" )

Figure 5

Radon soil-gas classification based on airborne radiometrics

  Col_Ramp <- RColorBrewer::brewer.pal(6, "Spectral")
  p1 <- ggplot() + 
    geom_sf(data = Ireland) +
    geom_sf(data = eU_Grids1km, aes(fill = Rn_class), colour = NA) +
    scale_fill_manual(name = "Rn Class",
                      values = Col_Ramp) +
    ggtitle("Radon soil-gas classification based on airborne radiometrics") 
 
  ggsave(p1,
         file = "Rresults/Figure_5.png",
         width = 14,
         height = 14,
         units = "cm" )

Figure 6

Radon Potential map based on gamma-ray spectrometry airborne radiometrics and subsoil permeability

  Col_Ramp <- RColorBrewer::brewer.pal(4, "Spectral")
  p1 <- ggplot() + 
    geom_sf(data = Ireland) +
    geom_sf(data = RP_Grids, aes(fill = RP_class), colour = NA) +
    scale_fill_manual(name = "GRP Class",
                      values = Col_Ramp) +
    ggtitle("Radon Potential map based on airborne radiometrics") 
  
   ggsave(p1,
          file = "Rresults/Figure_6.png",
         width = 14,
         height = 14,
         units = "cm" )

Figure 7

Histogram and probability plot of the new and old indoor radon measurements

  df_New <- as.data.frame(qqnorm(log10(InRn[InRn$Survey=="New",]$InRn),
                                 plot.it = FALSE))
  df_Old <- as.data.frame(qqnorm(log10(InRn[InRn$Survey=="Old",]$InRn),
                                 plot.it = FALSE))
   
  grDevices::png("Rresults/Figure_7.png",
                 width = 750,
                 height = 250,
                 pointsize = 16) 
  
    par(mfrow = c(1, 3))
    {hist(log10(InRn[InRn$Survey == "New",]$InRn),
           xlim = c(0.5, 4),
           ylim = c(0,1),
           freq = FALSE,
           col = 2, 
           main = "New survey",
           xlab = expression("log"[10] * " InRn" *  " [Bq" * m^-3 * " ]"))
     hist(log10(InRn[InRn$Survey=="Old",]$InRn),
           xlim = c(0.5, 4),
           ylim = c(0, 1),
           freq = FALSE,
           col = 2,
           main = "Old survey",
           xlab = expression("log"[10] * " InRn" *  " [Bq" * m^-3 * " ]"))
      
     plot(df_New$y,df_New$x,
           ylim = c(-4, 4),
           xlim = c(0,4), 
           axes = FALSE,
           frame.plot = TRUE,
           main = "Probability plot",
           xlab = "",
           ylab = "",
           pch = 19,
           col = "red",
           cex = 0.5)
      
     points(df_Old$y, df_Old$x, col = "black", pch = 19, cex = 0.4)
      
     ay<-c(0.01,0.1,1,5,20,50,80,95,99,99.9,99.999)
     at<-qnorm(ay/100)
     axis(2,labels=ay,at=at)
     mtext("Probability (%)", side= 2, line=3,cex=0.7)
      
     ax <- c(seq(0.1,1,0.1),
             seq(2,10,1),
             seq(20,100,10),
             seq(200,1000,100),
             seq(2000,10000,1000)
             )
     
     atx <- log10(ax)
     axis(1,labels = F, at = atx)
     
     ax <- 10^(seq(-4, 4, 1))
     atx <- log10(ax)
     axis(1, labels = ax, at = atx, col.ticks = 2)
      
     mtext(expression(" InRn" *  " [Bq" * m^-3 * " ]"),
           side = 1,
           line = 3,
           cex = 0.7)
     abline(v = atx, h = at, lty = 2, lwd = 0.1, col = "lightgray")
      
     legend("topleft",
            c("New survey","Old survey"),
            col = c("red","black"),
            pch = c(19,19),
            cex = 1,
            ncol = 1)
     }
  
  dev.off()

Figure 8

In-situ vs. airborne 226Ra concentration in soil (lines represent the trend lines of the percentiles 2.5th, 50th, 75th and 97.5th obtained in the simulation; the yellow points are the soil-gas radon samples; and the black points are the geometric mean of the radon measurements in each grid).

  p1 <- ggplot(eU_Grids1km) +
  geom_hline(yintercept = c(10, 30, 50, 70, 100),
             linetype = "dashed", color = "darkgrey", size = 0.5) +
  geom_point(aes(x = AM_Ra226, y = Rn_P025, colour = "purple"),
             size = 0.7,
             alpha = 0.1) +
  geom_point(aes(x = AM_Ra226, y = Rn_P50, colour = "green"),
             size = 0.7,
             alpha = 0.1) +
  geom_point(aes(x = AM_Ra226, y = Rn_P75, colour = "blue"),
             size = 0.7,
             alpha = 0.1) +
  geom_point(aes(x = AM_Ra226, y = Rn_P975, colour = "red"),
             size = 0.7,
             alpha = 0.1) +
  scale_colour_identity(name = "Percentiles",
                       breaks = c("purple", "green", "blue", "red"),
                       labels = c("2.5th", "50th", "75th", "97.5th"),
                       guide = "legend") +
  guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))) +
  geom_point(data = Soil_Rn, aes(x = AM_Ra226, y = SoilRn),
             size = 0.8,
             colour = "black") +
  geom_point(data = SG_Sum, aes(x = Ra_AM, y = Rn_GM),
             colour = "black",
             shape = 18,
             size = 5) + 
  xlim(0, 45) +
  ggtitle(
    "In-situ Rn-222 soil-gas vs. airborne Ra-226 (Grids 1km x 1km)"
  ) + 
  xlab(expression("Mean "*""^226*Ra * " [ " * "Bq "* kg^-1 * " ]")) + 
  ylab(expression(""^222*Rn * " [ " * "kBq "* m^-3 * " ]")) +
  scale_y_log10(breaks = breaks,
                minor_breaks = minor_breaks,
                limits = c(1, 1000)) 

  ggsave(p1,
         file = "Rresults/Figure_8.png",
         width = 20,
         height = 14,
         units = "cm" )

Figure 9

Boxplot of indoor radon measurements (InRn) in each Radon Potential class (L: Low; ML: Moderate-Low; MH: Moderate-High; and H: High). Boxes widths are proportional to the square-roots of the number of observations in each group

  p1 <- ggplot(InRn_RP, aes(x = RP_class, y = log10(InRn))) + 
    geom_boxplot(notch = TRUE, fill = "red") +
    ylab("log10(InRn)") +
    xlab("GRP")

  ggsave(p1,
         file = "Rresults/Figure_9.png",
         width = 7,
         height = 7,
         units = "cm" )